1 A well-defined HDX-MS experiment

This vignette describes how to analyse time-resolved differential HDX-MS experiments. The key elements are at least two conditions i.e. apo + antibody, apo + small molecule or protein closed + protien open, etc. The experiment can be replicated, though if there are sufficient time points analysed (>=3) then occasionally signficant results can be obtained. The data provided should be centroid-centric data. This package does not yet support analysis straight from raw spectra. Typically this will be provided as a .csv from tools such as dynamiX or HDExaminer.

2 Main elements of the package

The package relies of Bioconductor infrastructure so that it integrates with other data types and can benefit from advantages in other fields of mass-spectrometry. There are package specific object, classes and methods but importantly there is reuse of classes found in quantitative proteomics data, mainly the QFeatures object which extends the SummarisedExperiment class for mass spectrometry data. The focus of this package is on testing and visualisation of the testing results.

3 Data

We will begin with a structural variant experiment in which MHP and a structural variant were mixed in different proportions. HDX-MS was performed on these samples and we expect to see reproducible but subtle differences. We first load the data from the package and it is .csv format.

MBPpath <- system.file("extdata", "MBP.csv", package = "hdxstats")

We can now read in the .csv file and have a quick look at the .csv.

MBP <- read.csv(MBPpath)
head(MBP) # have a look
##   hx_sample pep_start pep_end pep_sequence pep_charge     d confidence  score
## 1       10%        19      30 VIWINGDKGYNG          2 2.120     medium 0.8686
## 2       10%        19      30 VIWINGDKGYNG          2 2.146     medium 0.8173
## 3       10%        19      30 VIWINGDKGYNG          2 2.143     medium 0.8839
## 4       10%        19      30 VIWINGDKGYNG          2    NA       <NA>     NA
## 5       10%        19      30 VIWINGDKGYNG          2    NA       <NA>     NA
## 6       10%        19      30 VIWINGDKGYNG          2    NA       <NA>     NA
##   hx_time time_unit replicate_cnt
## 1      30         s             1
## 2      30         s             2
## 3      30         s             3
## 4      30         s             4
## 5      30         s             5
## 6      30         s             6
length(unique(MBP$pep_sequence)) # peptide sequences
## [1] 115

Let us have a quick visualisation of some the data so that we can see some of the features

filter(MBP, pep_sequence == unique(MBP$pep_sequence[1]), pep_charge == 2) %>%
    ggplot(aes(x = hx_time, y = d, group = factor(replicate_cnt),
               color = factor(hx_sample,
                              unique(MBP$hx_sample)[c(7,5,1,2,3,4,6)]))) + 
    theme_classic() + geom_point(size = 2) + 
    scale_color_manual(values = brewer.pal(n = 7, name = "Set2")) + 
    labs(color = "experiment", x = "Deuterium Exposure", y = "Deuterium incoperation")
## Warning: Removed 96 rows containing missing values (geom_point).

We can see that the units of the time dimension are in seconds and that Deuterium incorporation has been normalized into Daltons.

4 Parsing to an object of class QFeatures

Working from a .csv is likely to cause issues downstream. Indeed, we run the risk of accidently changing the data or corrupting the file in some way. Secondly, all .csvs will be formatted slightly different and so making extensible tools for these files will be inefficient. Furthermore, working with a generic class used in other mass-spectrometry fields can speed up analysis and adoption of new methods. We will work the class QFeatures from the QFeatures class as it is a powerful and scalable way to store quantitative mass-spectrometry data.

Firstly, the data is storted in long format rather than wide format. We first switch the data to wide format.

MBP_wide <- pivot_wider(data.frame(MBP),
                        values_from = d,
                        names_from = c("hx_time", "replicate_cnt", "hx_sample"),
                        id_cols = c("pep_sequence", "pep_charge"))
head(MBP_wide)
## # A tibble: 6 × 198
##   pep_sequence   pep_c…¹ 30_1_…² 30_2_…³ 30_3_…⁴ 30_4_…⁵ 30_5_…⁶ 30_6_…⁷ 30_7_…⁸
##   <chr>            <int>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
## 1 VIWINGDKGYNG         2   2.12    2.15    2.14       NA      NA      NA      NA
## 2 VIWINGDKGYNGL        2   2.12    2.12    2.11       NA      NA      NA      NA
## 3 GDKGYNGLAEVG         3   0.552   0.555   0.553      NA      NA      NA      NA
## 4 LAEVGKKFEKDTG…       4   2.41    2.36    2.42       NA      NA      NA      NA
## 5 AEVGKKFEKDTGI…       4   0.458   0.425   0.573      NA      NA      NA      NA
## 6 TVEHPDKL             3   1.43    1.44    1.44       NA      NA      NA      NA
## # … with 189 more variables: `240_1_10%` <dbl>, `240_2_10%` <dbl>,
## #   `240_3_10%` <dbl>, `240_4_10%` <dbl>, `240_5_10%` <dbl>, `240_6_10%` <dbl>,
## #   `240_7_10%` <dbl>, `1800_1_10%` <dbl>, `1800_2_10%` <dbl>,
## #   `1800_3_10%` <dbl>, `1800_4_10%` <dbl>, `1800_5_10%` <dbl>,
## #   `1800_6_10%` <dbl>, `1800_7_10%` <dbl>, `14400_1_10%` <dbl>,
## #   `14400_2_10%` <dbl>, `14400_3_10%` <dbl>, `14400_4_10%` <dbl>,
## #   `14400_5_10%` <dbl>, `14400_6_10%` <dbl>, `14400_7_10%` <dbl>, …

We notice that there are many columns with NAs. The follow code chunk removes these columns.

MBP_wide <- MBP_wide[, colSums(is.na(MBP_wide)) != nrow(MBP_wide)]

We also note that the colnames are not very informative. We are going to format in a very specific way so that later functions can automatically infer the design from the column names. We provide in the format X(time)rep(replicate)cond(condition)

colnames(MBP_wide)[-c(1,2)]
##   [1] "30_1_10%"        "30_2_10%"        "30_3_10%"        "240_1_10%"      
##   [5] "240_2_10%"       "240_3_10%"       "1800_1_10%"      "1800_2_10%"     
##   [9] "1800_3_10%"      "14400_1_10%"     "14400_2_10%"     "14400_3_10%"    
##  [13] "30_1_15%"        "30_2_15%"        "30_3_15%"        "240_1_15%"      
##  [17] "240_2_15%"       "240_3_15%"       "1800_1_15%"      "1800_2_15%"     
##  [21] "1800_3_15%"      "14400_1_15%"     "14400_2_15%"     "14400_3_15%"    
##  [25] "30_1_20%"        "30_2_20%"        "30_3_20%"        "240_1_20%"      
##  [29] "240_2_20%"       "240_3_20%"       "1800_1_20%"      "1800_2_20%"     
##  [33] "1800_3_20%"      "14400_1_20%"     "14400_2_20%"     "14400_3_20%"    
##  [37] "30_1_25%"        "30_2_25%"        "30_3_25%"        "240_1_25%"      
##  [41] "240_2_25%"       "240_3_25%"       "1800_1_25%"      "1800_2_25%"     
##  [45] "1800_3_25%"      "14400_1_25%"     "14400_2_25%"     "14400_3_25%"    
##  [49] "30_1_5%"         "30_2_5%"         "30_3_5%"         "240_1_5%"       
##  [53] "240_2_5%"        "240_3_5%"        "1800_1_5%"       "1800_2_5%"      
##  [57] "1800_3_5%"       "14400_1_5%"      "14400_2_5%"      "14400_3_5%"     
##  [61] "30_1_W169G"      "30_2_W169G"      "30_3_W169G"      "240_1_W169G"    
##  [65] "240_2_W169G"     "240_3_W169G"     "1800_1_W169G"    "1800_2_W169G"   
##  [69] "1800_3_W169G"    "14400_1_W169G"   "14400_2_W169G"   "14400_3_W169G"  
##  [73] "30_1_WT Null"    "30_2_WT Null"    "30_3_WT Null"    "30_4_WT Null"   
##  [77] "30_5_WT Null"    "30_6_WT Null"    "30_7_WT Null"    "240_1_WT Null"  
##  [81] "240_2_WT Null"   "240_3_WT Null"   "240_4_WT Null"   "240_5_WT Null"  
##  [85] "240_6_WT Null"   "240_7_WT Null"   "1800_1_WT Null"  "1800_2_WT Null" 
##  [89] "1800_3_WT Null"  "1800_4_WT Null"  "1800_5_WT Null"  "1800_6_WT Null" 
##  [93] "1800_7_WT Null"  "14400_1_WT Null" "14400_2_WT Null" "14400_3_WT Null"
##  [97] "14400_4_WT Null" "14400_5_WT Null" "14400_6_WT Null" "14400_7_WT Null"
new.colnames <- gsub("0_", "0rep", paste0("X", colnames(MBP_wide)[-c(1,2)]))
new.colnames <- gsub("_", "cond", new.colnames)

# remove annoying % signs
new.colnames <- gsub("%", "", new.colnames)

# remove space (NULL could get confusing later and WT is clear)
new.colnames <- gsub(" .*", "", new.colnames)

new.colnames
##   [1] "X30rep1cond10"       "X30rep2cond10"       "X30rep3cond10"      
##   [4] "X240rep1cond10"      "X240rep2cond10"      "X240rep3cond10"     
##   [7] "X1800rep1cond10"     "X1800rep2cond10"     "X1800rep3cond10"    
##  [10] "X14400rep1cond10"    "X14400rep2cond10"    "X14400rep3cond10"   
##  [13] "X30rep1cond15"       "X30rep2cond15"       "X30rep3cond15"      
##  [16] "X240rep1cond15"      "X240rep2cond15"      "X240rep3cond15"     
##  [19] "X1800rep1cond15"     "X1800rep2cond15"     "X1800rep3cond15"    
##  [22] "X14400rep1cond15"    "X14400rep2cond15"    "X14400rep3cond15"   
##  [25] "X30rep1cond20"       "X30rep2cond20"       "X30rep3cond20"      
##  [28] "X240rep1cond20"      "X240rep2cond20"      "X240rep3cond20"     
##  [31] "X1800rep1cond20"     "X1800rep2cond20"     "X1800rep3cond20"    
##  [34] "X14400rep1cond20"    "X14400rep2cond20"    "X14400rep3cond20"   
##  [37] "X30rep1cond25"       "X30rep2cond25"       "X30rep3cond25"      
##  [40] "X240rep1cond25"      "X240rep2cond25"      "X240rep3cond25"     
##  [43] "X1800rep1cond25"     "X1800rep2cond25"     "X1800rep3cond25"    
##  [46] "X14400rep1cond25"    "X14400rep2cond25"    "X14400rep3cond25"   
##  [49] "X30rep1cond5"        "X30rep2cond5"        "X30rep3cond5"       
##  [52] "X240rep1cond5"       "X240rep2cond5"       "X240rep3cond5"      
##  [55] "X1800rep1cond5"      "X1800rep2cond5"      "X1800rep3cond5"     
##  [58] "X14400rep1cond5"     "X14400rep2cond5"     "X14400rep3cond5"    
##  [61] "X30rep1condW169G"    "X30rep2condW169G"    "X30rep3condW169G"   
##  [64] "X240rep1condW169G"   "X240rep2condW169G"   "X240rep3condW169G"  
##  [67] "X1800rep1condW169G"  "X1800rep2condW169G"  "X1800rep3condW169G" 
##  [70] "X14400rep1condW169G" "X14400rep2condW169G" "X14400rep3condW169G"
##  [73] "X30rep1condWT"       "X30rep2condWT"       "X30rep3condWT"      
##  [76] "X30rep4condWT"       "X30rep5condWT"       "X30rep6condWT"      
##  [79] "X30rep7condWT"       "X240rep1condWT"      "X240rep2condWT"     
##  [82] "X240rep3condWT"      "X240rep4condWT"      "X240rep5condWT"     
##  [85] "X240rep6condWT"      "X240rep7condWT"      "X1800rep1condWT"    
##  [88] "X1800rep2condWT"     "X1800rep3condWT"     "X1800rep4condWT"    
##  [91] "X1800rep5condWT"     "X1800rep6condWT"     "X1800rep7condWT"    
##  [94] "X14400rep1condWT"    "X14400rep2condWT"    "X14400rep3condWT"   
##  [97] "X14400rep4condWT"    "X14400rep5condWT"    "X14400rep6condWT"   
## [100] "X14400rep7condWT"

We will now parse the data into an object of class QFeatures, we have provided a function to assist with this in the package. If you want to do this yourself use the readQFeatures function from the QFeatures package.

MBPqDF <- parseDeutData(object = DataFrame(MBP_wide),
                        design = new.colnames,
                        quantcol = 3:102)

5 Heatmap visualisations of HDX data

To help us get used to the QFeatures we show how to generate a heatmap of these data from this object:

pheatmap(t(assay(MBPqDF)),
         cluster_rows = FALSE, 
         cluster_cols = FALSE,
         color = brewer.pal(n = 9, name = "BuPu"),
         main = "Stuctural variant deuterium incoperation heatmap", 
         fontsize = 14,
         legend_breaks = c(0, 2, 4, 6, 8, 10, 12, max(assay(MBPqDF))),
         legend_labels = c("0", "2", "4", "6", "8","10", "12", "Incorporation"))

If you prefer to have the start-to-end residue numbers in the heatmap instead you can change the plot as follows:

regions <- unique(MBP[,c("pep_start", "pep_end")])
xannot <- paste0("[", regions[,1], ",", regions[,2], "]")
pheatmap(t(assay(MBPqDF)),
         cluster_rows = FALSE, 
         cluster_cols = FALSE,
         color = brewer.pal(n = 9, name = "BuPu"),
         main = "Stuctural variant deuterium incoperation heatmap", 
         fontsize = 14,
         legend_breaks = c(0, 2, 4, 6, 8, 10, 12, max(assay(MBPqDF))),
         legend_labels = c("0", "2", "4", "6", "8","10", "12", "Incorporation"),
         labels_col = xannot)

It maybe useful to normalize HDX-MS data for either interpretation or visualization purposes. We can normalize by the number of exchangeable amides or by using back-exchange correction values. We first use percentage incorporation as normalisation and visualise as a heatmap.

MBPqDF_norm1 <- normalisehdx(MBPqDF,
                             sequence = unique(MBP$pep_sequence),
                             method = "pc")


pheatmap(t(assay(MBPqDF_norm1)),
         cluster_rows = FALSE, 
         cluster_cols = FALSE,
         color = brewer.pal(n = 9, name = "BuPu"),
         main = "Stuctural variant deuterium incoperation heatmap normalised", 
         fontsize = 14,
         legend_breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1, 1.2),
         legend_labels = c("0", "0.2", "0.4", "0.6", "0.8","1", "Incorporation"),
         labels_col = xannot)

Now, we demonstrate a back-exchange correction calculation. The back-exchange value are fictious by the code chunk below demonstrates how to set this up.

# made-up correction factor
correction <- (exchangeableAmides(unique(MBP$pep_sequence)) + 1) * 0.9

MBPqDF_norm2 <- normalisehdx(MBPqDF,
                             sequence = unique(MBP$pep_sequence),
                             method = "bc", 
                             correction = correction)


pheatmap(t(assay(MBPqDF_norm2)),
         cluster_rows = FALSE, 
         cluster_cols = FALSE,
         color = brewer.pal(n = 9, name = "BuPu"),
         main = "Stuctural variant deuterium incoperation heatmap normalised", 
         fontsize = 14,
         legend_breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1, 1.2),
         legend_labels = c("0", "0.2", "0.4", "0.6", "0.8","1", "Incorporation"),
         labels_col = xannot)

6 Functional data analysis of HDX-MS data

The hdxstats package uses an empirical Bayes functional approach to analyse the data. We explain this idea in steps so that we can get an idea of the approach. First we fit the parametric model to the data. This will allow us to explore the HdxStatModel class.

res <- differentialUptakeKinetics(object = MBPqDF[,1:100], #provide a QFeature object
                                  feature = rownames(MBPqDF)[[1]][37], # which peptide to do we fit
                                  start = list(a = NULL, b = NULL,  d = NULL, p = 1)) # what are the starting parameter guesses
## Warning in differentialUptakeKinetics(object = MBPqDF[, 1:100], feature =
## rownames(MBPqDF)[[1]][37], : NAs introduced by coercion

Here, we see the HdxStatModel class, and that a Functional Model was applied to the data and a total of 7 models were fitted.

res
## Object of class "HdxStatModel"
## Method: Functional Model 
## Fitted 7

The nullmodel and alternative slots of an instance of HdxStatModel provide the underlying fitted models. The method and formula slots provide vital information about what analysis was performed. The vis slot provides a ggplot object so that we can visualise the functional fits.

res@vis

Since this is a ggplot object, we can customise in the usual grammatical ways.

res@vis + scale_color_manual(values = brewer.pal(n = 8, name = "Set2"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

A number of standard methods are available and can be applied to a HdxStatModels, these extend the usual base stats methods. These include

  1. anova: An analysis of variance
  2. logLik: The log-likelihood of all the fitted models
  3. residuals: The residuals for the fitted models
  4. vcov: The variance-covariance matrix between parameters of the models
  5. likRatio: The likelihood ratio between null and alternative models
  6. wilk: Applies wilk’s theorem to obtain a p-value from the liklihood ratio
  7. coef: The fitted model coefficients
  8. deviance: The deviance of the fitted models
  9. summary: The statistical summary of the models.
anova(res)
## Analysis of Variance Table
## 
## Model 1: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
## Model 2: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
## Model 3: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
## Model 4: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
## Model 5: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
## Model 6: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
## Model 7: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
## Model 8: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
##   Res.Df Res.Sum Sq  Df  Sum Sq F value    Pr(>F)    
## 1     96    14.0842                                  
## 2      8     0.0788  88 14.0054 16.1547 0.0001443 ***
## 3      8     0.0969   0  0.0000                      
## 4      8     0.0622   0  0.0000                      
## 5      8     0.0439   0  0.0000                      
## 6      8     0.0891   0  0.0000                      
## 7      8     0.1381   0  0.0000                      
## 8     24     0.1978 -16 -0.0597  0.2162 0.9955201    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
logLik(res)
##       null       alt1       alt2       alt3       alt4       alt5       alt6 
## -43.887971  13.126183  11.888973  14.550061  16.637751  12.387702   9.762426 
##       alt7 
##  29.610172
residuals(res)
## $null
##   [1] -0.104877060 -0.200877060 -0.116877060 -0.041071747 -0.067071747
##   [6] -0.005071747 -0.261216241 -0.227216241 -0.299216241 -0.048625259
##  [11] -0.133625259 -0.110625259 -0.066877060 -0.093877060 -0.179877060
##  [16] -0.007071747  0.025928253  0.025928253 -0.271216241 -0.155216241
##  [21] -0.190216241  0.032374741 -0.015625259  0.059374741 -0.011877060
##  [26] -0.101877060 -0.010877060  0.036928253  0.022928253  0.054928253
##  [31] -0.167216241 -0.130216241 -0.184216241  0.002374741  0.065374741
##  [36] -0.001625259 -0.020877060 -0.071877060 -0.135877060  0.006928253
##  [41]  0.038928253  0.061928253 -0.081216241 -0.097216241 -0.080216241
##  [46] -0.065625259 -0.014625259  0.034374741 -0.109877060 -0.097877060
##  [51] -0.100877060 -0.032071747 -0.064071747 -0.020071747 -0.219216241
##  [56] -0.269216241 -0.371216241 -0.086625259  0.004374741 -0.067625259
##  [61]  0.394122940  0.451122940  0.430122940  1.206928253  1.253928253
##  [66]  1.193928253  1.312783759  1.230783759  1.314783759  0.676374741
##  [71]  0.673374741  0.659374741 -0.164877060 -0.120877060 -0.115877060
##  [76] -0.119877060 -0.106877060 -0.231877060 -0.175877060 -0.114071747
##  [81] -0.196071747 -0.144071747 -0.137071747 -0.154071747 -0.119071747
##  [86] -0.192071747 -0.467216241 -0.448216241 -0.378216241 -0.401216241
##  [91] -0.476216241 -0.416216241 -0.419216241 -0.216625259 -0.197625259
##  [96] -0.128625259 -0.108625259 -0.137625259 -0.192625259 -0.097625259
## attr(,"label")
## [1] "Residuals"
## 
## $alt1
##  [1] -0.009899831 -0.105899831 -0.021899831  0.106886146  0.080886146
##  [6]  0.142886146 -0.081792295 -0.047792295 -0.119792295  0.070965380
## [11] -0.014034620  0.008965380
## attr(,"label")
## [1] "Residuals"
## 
## $alt2
##  [1] -0.004344068 -0.031344068 -0.117344068  0.099861189  0.132861189
##  [6]  0.132861189 -0.157170607 -0.041170607 -0.076170607  0.031008519
## [11] -0.016991481  0.058008519
## attr(,"label")
## [1] "Residuals"
## 
## $alt3
##  [1] -1.156564e-02 -1.015656e-01 -1.056564e-02  9.740613e-02  8.340613e-02
##  [6]  1.154061e-01 -8.042102e-02 -4.342102e-02 -9.742102e-02  9.344649e-05
## [11]  6.309345e-02 -3.906554e-03
## attr(,"label")
## [1] "Residuals"
## 
## $alt4
##  [1]  0.01874941 -0.03225059 -0.09625059  0.04059671  0.07259671  0.09559671
##  [7] -0.05577593 -0.07177593 -0.05477593 -0.03282707  0.01817293  0.06717293
## attr(,"label")
## [1] "Residuals"
## 
## $alt5
##  [1] -0.056594876 -0.044594876 -0.047594876  0.123090620  0.091090620
##  [6]  0.135090620 -0.018095673 -0.068095673 -0.170095673 -0.014489369
## [11]  0.076510631  0.004510631
## attr(,"label")
## [1] "Residuals"
## 
## $alt6
##  [1] -0.11280750 -0.05580750 -0.07680750  0.14555184  0.19255184  0.13255184
##  [7] -0.07978353 -0.16178353 -0.07778353  0.03073823  0.02773823  0.01373823
## attr(,"label")
## [1] "Residuals"
## 
## $alt7
##  [1] -0.066801632 -0.022801632 -0.017801632 -0.021801632 -0.008801632
##  [6] -0.133801632 -0.077801632  0.150954795  0.068954795  0.120954795
## [11]  0.127954795  0.110954795  0.145954795  0.072954795 -0.119154441
## [16] -0.100154441 -0.030154441 -0.053154441 -0.128154441 -0.068154441
## [21] -0.071154441 -0.042026225 -0.023026225  0.045973775  0.065973775
## [26]  0.036973775 -0.018026225  0.076973775
## attr(,"label")
## [1] "Residuals"
vcov(res)
## $nullvcov
##             a            b            d            p
## a 63688.43320 -46.23397301 -626.8125656 -62.21085974
## b   -46.23397   0.03365449    0.4503213   0.04493626
## d  -626.81257   0.45032134    6.4159444   0.62375397
## p   -62.21086   0.04493626    0.6237540   0.06132549
## 
## $altvcov
## $altvcov[[1]]
##              a             b             d             p
## a 6742815.2868 -547.82144203 -4769.5078374 -537.38725418
## b    -547.8214    0.04451203     0.3867961    0.04362297
## d   -4769.5078    0.38679609     3.4950552    0.38635258
## p    -537.3873    0.04362297     0.3863526    0.04316242
## 
## $altvcov[[2]]
##              a             b             d             p
## a 7313985.8084 -583.88942708 -5225.6175495 -601.78890198
## b    -583.8894    0.04661769     0.4163713    0.04799888
## d   -5225.6175    0.41637133     3.8728898    0.43724610
## p    -601.7889    0.04799888     0.4372461    0.04991269
## 
## $altvcov[[3]]
##              a             b             d             p
## a 2578698.4241 -265.74265461 -2416.1108601 -286.02297154
## b    -265.7427    0.02738946     0.2484133    0.02944372
## d   -2416.1109    0.24841326     2.3505629    0.27265279
## p    -286.0230    0.02944372     0.2726528    0.03198702
## 
## $altvcov[[4]]
##             a            b            d           p
## a 5878.837748 -5.771596667 -98.50742036 -9.91292316
## b   -5.771597  0.005707975   0.09499339  0.00964856
## d  -98.507420  0.094993390   1.72282460  0.16951115
## p   -9.912923  0.009648560   0.16951115  0.01688395
## 
## $altvcov[[5]]
##              a             b             d             p
## a 7678516.7393 -473.78322396 -4244.5997729 -579.41593289
## b    -473.7832    0.02923634     0.2613853    0.03571877
## d   -4244.5998    0.26138533     2.4457797    0.32634826
## p    -579.4159    0.03571877     0.3263483    0.04411006
## 
## $altvcov[[6]]
##            a            b          d            p
## a  3.0294730  0.137107637 -2.1842144 -0.144502501
## b  0.1371076  0.006320457 -0.1003236 -0.006582307
## d -2.1842144 -0.100323607  1.6073619  0.103919672
## p -0.1445025 -0.006582307  0.1039197  0.006960797
## 
## $altvcov[[7]]
##              a             b             d             p
## a 2659633.5070 -1.196868e+02 -1.101485e+03 -1.857941e+02
## b    -119.6868  5.386547e-03  4.946382e-02  8.352965e-03
## d   -1101.4847  4.946382e-02  4.789887e-01  7.861826e-02
## p    -185.7941  8.352965e-03  7.861826e-02  1.310929e-02
likRatio(res)
##    logLR 
## 303.7025
wilk(res)
##      p-value 
## 3.011008e-50
coef(res)
##               a           b          d         p
## null  33.440379 0.031873046 0.00000000 0.2088560
## alt1 105.742817 0.009119417 0.06708871 0.2079068
## alt2 105.112047 0.008921838 0.12760063 0.2117928
## alt3  90.536895 0.010024700 0.23941111 0.2140327
## alt4  26.480450 0.038936619 0.00000000 0.2153552
## alt5 108.655735 0.007103991 0.39105083 0.2265999
## alt6   8.328166 0.131389103 0.00000000 0.3086935
## alt7 112.222583 0.005336975 0.62105325 0.2473897
deviance(res)
##        null        alt1        alt2        alt3        alt4        alt5 
## 14.08418504  0.07881376  0.09686221  0.06216392  0.04389614  0.08913641 
##        alt6        alt7 
##  0.13806353  0.19776322
summary(res)
## $null
## 
## Formula: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)
## a  33.44038  252.36567   0.133    0.895
## b   0.03187    0.18345   0.174    0.862
## d   0.00000    2.53297   0.000    1.000
## p   0.20886    0.24764   0.843    0.401
## 
## Residual standard error: 0.383 on 96 degrees of freedom
## 
## Number of iterations to convergence: 45 
## Achieved convergence tolerance: 1.49e-08
## 
## 
## $alt1
## 
## Formula: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)
## a 1.057e+02  2.597e+03   0.041    0.969
## b 9.119e-03  2.110e-01   0.043    0.967
## d 6.709e-02  1.870e+00   0.036    0.972
## p 2.079e-01  2.078e-01   1.001    0.346
## 
## Residual standard error: 0.09926 on 8 degrees of freedom
## 
## Number of iterations till stop: 96 
## Achieved convergence tolerance: 1.49e-08
## Reason stopped: Number of calls to `fcn' has reached or exceeded `maxfev' == 500.
## 
## 
## $alt2
## 
## Formula: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)
## a 1.051e+02  2.704e+03   0.039    0.970
## b 8.922e-03  2.159e-01   0.041    0.968
## d 1.276e-01  1.968e+00   0.065    0.950
## p 2.118e-01  2.234e-01   0.948    0.371
## 
## Residual standard error: 0.11 on 8 degrees of freedom
## 
## Number of iterations till stop: 95 
## Achieved convergence tolerance: 1.49e-08
## Reason stopped: Number of calls to `fcn' has reached or exceeded `maxfev' == 500.
## 
## 
## $alt3
## 
## Formula: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)
## a 9.054e+01  1.606e+03   0.056    0.956
## b 1.002e-02  1.655e-01   0.061    0.953
## d 2.394e-01  1.533e+00   0.156    0.880
## p 2.140e-01  1.789e-01   1.197    0.266
## 
## Residual standard error: 0.08815 on 8 degrees of freedom
## 
## Number of iterations till stop: 95 
## Achieved convergence tolerance: 1.49e-08
## Reason stopped: Number of calls to `fcn' has reached or exceeded `maxfev' == 500.
## 
## 
## $alt4
## 
## Formula: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)
## a 26.48045   76.67358   0.345    0.739
## b  0.03894    0.07555   0.515    0.620
## d  0.00000    1.31256   0.000    1.000
## p  0.21536    0.12994   1.657    0.136
## 
## Residual standard error: 0.07407 on 8 degrees of freedom
## 
## Number of iterations to convergence: 45 
## Achieved convergence tolerance: 1.49e-08
## 
## 
## $alt5
## 
## Formula: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)
## a 1.087e+02  2.771e+03   0.039    0.970
## b 7.104e-03  1.710e-01   0.042    0.968
## d 3.911e-01  1.564e+00   0.250    0.809
## p 2.266e-01  2.100e-01   1.079    0.312
## 
## Residual standard error: 0.1056 on 8 degrees of freedom
## 
## Number of iterations till stop: 95 
## Achieved convergence tolerance: 1.49e-08
## Reason stopped: Number of calls to `fcn' has reached or exceeded `maxfev' == 500.
## 
## 
## $alt6
## 
## Formula: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)   
## a  8.32817    1.74054   4.785  0.00138 **
## b  0.13139    0.07950   1.653  0.13700   
## d  0.00000    1.26782   0.000  1.00000   
## p  0.30869    0.08343   3.700  0.00604 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1314 on 8 degrees of freedom
## 
## Number of iterations to convergence: 32 
## Achieved convergence tolerance: 1.49e-08
## 
## 
## $alt7
## 
## Formula: value ~ a * (1 - exp(-b * (timepoint)^p)) + d
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)  
## a 1.122e+02  1.631e+03   0.069   0.9457  
## b 5.337e-03  7.339e-02   0.073   0.9426  
## d 6.211e-01  6.921e-01   0.897   0.3784  
## p 2.474e-01  1.145e-01   2.161   0.0409 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09078 on 24 degrees of freedom
## 
## Number of iterations till stop: 95 
## Achieved convergence tolerance: 1.49e-08
## Reason stopped: Number of calls to `fcn' has reached or exceeded `maxfev' == 500.

7 Analysis of a typical HDX-MS experiment

We have seen the basic aspects of our functional modelling approach. We now wish to roll out our method across all peptides in the experiment. The fitUptakeKinetics function allows us to apply our modelling approach across all the peptide in the experiment. We need to provide a QFeatures object and the features for which we are fitting the model. The design will be extracted from the column names or you can provide a design yourself. The parameter initilisation should also be provided. Sometimes the model can’t be fit on the kinetics. This is either because there is not enough data or through lack of convergence. An error will be reported in these cases but this should not perturb the user. You may wish to try a few starting values if there excessive models that fail fitting.

res <- fitUptakeKinetics(object = MBPqDF[,c(1:24)],
                         feature = rownames(MBPqDF[,c(1:24)])[[1]],
                         start = list(a = NULL, b = 0.001,  d = NULL, p = 1))
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"

The code chunk above returns a class HdxStatModels indicating that a number of models for peptide have been fit. This is simply a holder for a list of HdxStatModel instances.

res
## Object of class "HdxStatModels"
## Number of models 102

We can easily examine indivual fits by going to the underyling HdxStatModel class:

res@statmodels[[1]]@vis + scale_color_manual(values = brewer.pal(n = 2, name = "Set2"))
## Warning in brewer.pal(n = 2, name = "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

We now wish to apply statistical analysis to these fitted curves. Our approach is an empirical Bayes testing procedure, which borrows information across peptides to stablise variance estimates. Here, we need to provide the original data that was analysed and the HdxStatModels class. The following code chunk returns an object of class HdxStatRes. This object tell us that statistical analysis was performed using our Functional model.

out <- processFunctional(object = MBPqDF[,1:24], params = res)
## Warning in stats::optim(x = c(0.00109178357660077, 4.02437669894543e-05, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in densfun(x, parm, ...): NaNs produced

## Warning in densfun(x, parm, ...): NaNs produced

## Warning in densfun(x, parm, ...): NaNs produced
## Warning in stats::optim(x = c(0.00153823834115607, 0.00188329333146131, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in densfun(x, parm, ...): NaNs produced
out
## Object of class "HdxStatRes"
## Analysed using Functional model

The main slot of interest is the results slot which returns quantities of interest such as p-values and fdr corrected p-values because of multiple testing. The following is the DataFrame of interest.

out@results
## DataFrame with 102 rows and 9 columns
##                           Fstat.Fstat Fstat.numerator Fstat.denomenator
##                                <list>          <list>            <list>
## VIWINGDKGYNG_2               0.709762      0.00109178        0.00153824
## VIWINGDKGYNGL_2             0.0213688     4.02438e-05        0.00188329
## GDKGYNGLAEVG_3               0.980103      0.00129713        0.00132346
## LAEVGKKFEKDTGIKVTVEHPDK_4    0.809885      0.00117025        0.00144496
## AEVGKKFEKDTGIKV_4            0.806449      0.00140588        0.00174329
## ...                               ...             ...               ...
## YAVRTAVINAASGRQTVDE_3        -3.96761        -2.79196          0.703687
## AVINA_1                      -3.87726       -0.213395         0.0550377
## AVINAASGRQTVDEAL_2           0.905874      0.00277575        0.00306417
## ASGRQTVDE_2                   0.37392      0.00288595        0.00771809
## ASGRQTVDEAL_2                0.364883      0.00215214        0.00589818
##                               pvals       fdr ebayes.pvals ebayes.fdr
##                           <numeric> <numeric>    <numeric>  <numeric>
## VIWINGDKGYNG_2             0.597056         1     0.589127          1
## VIWINGDKGYNGL_2            0.999008         1     0.998946          1
## GDKGYNGLAEVG_3             0.445918         1     0.442886          1
## LAEVGKKFEKDTGIKVTVEHPDK_4  0.536935         1     0.530748          1
## AEVGKKFEKDTGIKV_4          0.538920         1     0.526604          1
## ...                             ...       ...          ...        ...
## YAVRTAVINAASGRQTVDE_3      1.000000         1     1.000000          1
## AVINA_1                    1.000000         1     1.000000          1
## AVINAASGRQTVDEAL_2         0.483850         1     0.457787          1
## ASGRQTVDE_2                0.823819         1     0.803136          1
## ASGRQTVDEAL_2              0.830029         1     0.811177          1
##                           empirical.fdr fitcomplete
##                               <numeric>   <integer>
## VIWINGDKGYNG_2                 0.727165           1
## VIWINGDKGYNGL_2                0.762892           2
## GDKGYNGLAEVG_3                 0.727165           3
## LAEVGKKFEKDTGIKVTVEHPDK_4      0.727165           4
## AEVGKKFEKDTGIKV_4              0.727165           5
## ...                                 ...         ...
## YAVRTAVINAASGRQTVDE_3          1.000000          98
## AVINA_1                        1.000000          99
## AVINAASGRQTVDEAL_2             0.727165         100
## ASGRQTVDE_2                    0.727165         101
## ASGRQTVDEAL_2                  0.727165         102

We can now examine the peptides for which the false discovery rate is less than 0.05

which(out@results$ebayes.fdr < 0.05)
## EAAFNKGETAM_2   AAFNKGETA_2  AAFNKGETAM_2 WYAVRTAVINA_2 
##            55            56            57            96

Let us visualise some of these examples:

res@statmodels[[42]]@vis + res@statmodels[[45]]@vis

As we can see our model has picked up some subtle differences, we can further visualise these using a forest plot. We can see the the functions are very similar as the parameters are almost identical (a,b,p,d). However, we can see that the deuterium differences are lower in 10% structural variant condition.

fp <- forestPlot(params = res@statmodels[[42]])

We can produce a table to actual numbers. We see that at all 4 timepoints the deuterium difference is negative, though the confidence intervals overlap with 0. Our functional approach is picking up this small but reproducible difference.

knitr::kable(fp$data)
Estimate confL confU rownames condition
a 11.1548050 -57.0931637 79.4027737 a 1
b 1.2581965 -4.0590049 6.5753979 b 1
d 0.0000000 -67.8160474 67.8160474 d 1
p 0.1332137 -0.2441360 0.5105634 p 1
a1 11.1204905 -61.3165379 83.5575188 a 2
b1 1.3563310 -4.3737136 7.0863756 b 2
d1 0.0000000 -72.0481065 72.0481065 d 2
p1 0.1268633 -0.2395994 0.4933261 p 2
1 -0.1336667 -0.2728705 0.0055372 Timepoint 30 Deuterium Difference
2 -0.0350000 -0.2080751 0.1380751 Timepoint 240 Deuterium Difference
3 -0.0196667 -0.0873897 0.0480564 Timepoint 1800 Deuterium Difference
4 0.0343333 -0.1120953 0.1807620 Timepoint 14400 Deuterium Difference

It is also possible to visualize, these plots on a different scale. Of course, changing the natural scaling will emphasis different parts of the plot and could distort interpretation. In particular, if a log transform is used then care should be taken when interpreting values around 0. We suggest examining the numerical values in a forest plot or table alongside any transformation of the variables. We suggest using the pseudo log transform as this allows control the linearity of the plot, clearly demonstrating this a choice of visualisation (and not of statistical modelling). The parameter sigma below controls the scaling factor of the linear part of the transformation.

res@statmodels[[42]]@vis + scale_x_continuous(
    trans = pseudo_log_trans(base = 10, sigma = 0.01), breaks = c(0, 10^(1:7)))

res@statmodels[[42]]@vis + scale_x_continuous(
    trans = pseudo_log_trans(base = 10, sigma = 0.0001), breaks = c(0, 10^(1:7)))

res@statmodels[[42]]@vis + scale_x_continuous(
    trans = pseudo_log_trans(base = 10, sigma = 10), breaks = c(0, 10^(1:7)))

Let’s us now have a look a situation where the changes are more dramatic.

res_wt <- fitUptakeKinetics(object = MBPqDF[, c(61:100)],
                            feature = rownames(MBPqDF[, c(61:100)])[[1]],
                            start = list(a = NULL, b = 0.001,  d = NULL, p = 1))
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion

## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
out_wt <- processFunctional(object = MBPqDF[, c(61:100)], params = res_wt)
## Warning in stats::optim(x = c(0.0840870774115754, 0.0705385629222093, 0.0354279988095238, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in stats::optim(x = c(0.00199791426867962, 0.00247716448153206, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in densfun(x, parm, ...): NaNs produced

We can visualise some of the result and generate plots.

res_wt@statmodels[[27]]@vis/res_wt@statmodels[[28]]@vis + plot_layout(guides = "collect")|(forestPlot(params = res_wt@statmodels[[27]], condition = c("WT", "W169G"))/forestPlot(params = res_wt@statmodels[[28]], condition = c("WT", "W169G")) + plot_layout(guides = "collect")) + 
    plot_annotation(tag_levels = 'a') +  plot_layout(widths = c(1, 1))

8 An epitope mapping experiment

We now describe the analysis of an epitope mapping experiment. Here, the data analysis is more challenging, since only 1 replicate in each condition, apo and antibody, was performed. If we make some simplifying assumptions rigorous statistical analysis can still be performed.

The experiment was performed on HOIP-RBR, we loaded the data below from inside the package

HOIPpath <- system.file("extdata", "N64184_1a2_state.csv", package = "hdxstats")
HOIP <- read.csv(HOIPpath)
unique(HOIP$State)
##  [1] "apo"     "dAb13_1" "dAb13_2" "dAb25_1" "dAb25_2" "dAb27_1" "dAb27_2"
##  [8] "dAb2_1"  "dAb2_2"  "dAb6_1"  "dAb6_2"
HOIP$Exposure <- HOIP$Exposure * 60 #convert to seconds
filter(HOIP, Sequence == unique(HOIP$Sequence[1])) %>%
    ggplot(aes(x = Exposure,
               y = Center,
               color = factor(State, unique(HOIP$State)))) +
    theme_classic() + geom_point(size = 3) + 
    scale_color_manual(values = colorRampPalette(brewer.pal(8, name = "Set2"))(11)) + 
    labs(color = "experiment", x = "Deuterium Exposure", y = "Deuterium incoperation")

As before we need to convert data to an object of classes QFeatures for ease of analysis.

First, we put the data into a DataFrame object. Currently, its in long format so we switch to a wide format

HOIP_wide <- pivot_wider(data.frame(HOIP),
                         values_from = Center,
                         names_from = c("Exposure", "State"),
                         id_cols = c("Sequence"))

Now remove all columns with only NAs

HOIP_wide <- HOIP_wide[, colSums(is.na(HOIP_wide)) != nrow(HOIP_wide)]

The colanmes are not very informative, provide in the format X(time)rep(repliate)cond(condition)

colnames(HOIP_wide)[-c(1)]
##  [1] "0_apo"       "30_apo"      "300_apo"     "0_dAb13_1"   "30_dAb13_1" 
##  [6] "300_dAb13_1" "0_dAb13_2"   "30_dAb13_2"  "300_dAb13_2" "0_dAb25_1"  
## [11] "30_dAb25_1"  "300_dAb25_1" "0_dAb25_2"   "30_dAb25_2"  "300_dAb25_2"
## [16] "0_dAb27_1"   "30_dAb27_1"  "300_dAb27_1" "0_dAb27_2"   "30_dAb27_2" 
## [21] "300_dAb27_2" "0_dAb2_1"    "30_dAb2_1"   "300_dAb2_1"  "0_dAb2_2"   
## [26] "30_dAb2_2"   "300_dAb2_2"  "0_dAb6_1"    "30_dAb6_1"   "300_dAb6_1" 
## [31] "0_dAb6_2"    "30_dAb6_2"   "300_dAb6_2"
new.colnames <- gsub("0_", "0rep1", paste0("X", colnames(HOIP_wide)[-c(1)]))
new.colnames <- gsub("rep1", "rep1cond", new.colnames)

# remove annoying % signs
new.colnames <- gsub("%", "", new.colnames)

# remove space (NULL could get confusing later and WT is clear)
new.colnames <- gsub(" .*", "", new.colnames)

Now, we can provide rownames and convert the data to a QFeatures object:

qDF <- parseDeutData(object = DataFrame(HOIP_wide),
                     design = new.colnames,
                     quantcol = 2:34,
                     rownames = HOIP_wide$Sequence)

As before, we can produce a heatmap, we perform a simple normalisation for ease of visualisation:

mat <- assay(qDF)
mat <- apply(mat, 2, function(x) x - assay(qDF)[,1])

pheatmap(t(mat),
         cluster_rows = FALSE,
         cluster_cols = FALSE,
         color = brewer.pal(n = 9, name = "BuPu"),
         main = "HOIP RBR heatmap",
         fontsize = 14,
         legend_breaks = c(0, 2, 4, 6,8,10,12, max(assay(qDF))),
         legend_labels = c("0", "2", "4", "6", "8","10", "12", "Incorporation"))

Let us first perform a quick test:

res <- differentialUptakeKinetics(object = qDF[,1:33],
                                  feature = rownames(qDF)[[1]][3],
                                  start = list(a = NULL, b = 0.01,  d = NULL),
                                  formula = value ~ a * (1 - exp(-b*(timepoint))) + d) 
## Warning in differentialUptakeKinetics(object = qDF[, 1:33], feature =
## rownames(qDF)[[1]][3], : NAs introduced by coercion
res@vis+ scale_color_manual(values = colorRampPalette(brewer.pal(8, name = "Set2"))(11)) 
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

Whilst this analysis performs good fits for the functions, there are too many degrees of freedom to perform sound statistical analysis. Hence, we normalize to remove the degree of freedom for the intercept. For simplicity and to preserve the original matrix, we reprocess the data. We then fit a simplified kinetic model, where only the plateau is inferred.

cn <- new.colnames[c(1:3,10:12)]
HOIP_wide_nrm <- data.frame(HOIP_wide)
HOIP_wide_nrm[, c(2:4)] <- HOIP_wide_nrm[,c(2:4)] - HOIP_wide_nrm[,c(2)] # normalise by intercept
HOIP_wide_nrm[, c(11:13)] <- HOIP_wide_nrm[,c(11:13)] - HOIP_wide_nrm[,c(11)] # normalised by intercept

newqDF <- parseDeutData(object = DataFrame(HOIP_wide_nrm),
                        design = cn,
                        quantcol = c(2:4, 11:13), rownames = HOIP_wide$Sequence)

res_all <- fitUptakeKinetics(object = newqDF[,1:6],
                             feature = rownames(newqDF[,1:6])[[1]],
                             start = list(a = NULL),
                             formula = value ~ a * (1 - exp(-0.05*(timepoint))), maxAttempts = 1)
  
funresdAb25_1 <- processFunctional(object = newqDF[,1:6],
                                   params = res_all)
## Warning in stats::optim(x = c(0.0236132769138663, 0.0111598488727292, 0.130743883316999, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in stats::optim(x = c(0.0466717805074177, 0.0136102073076068, 0.041005949389079, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in densfun(x, parm, ...): NaNs produced

We can have a look at the results:

funresdAb25_1@results
## DataFrame with 110 rows and 9 columns
##                           Fstat.Fstat Fstat.numerator Fstat.denomenator
##                                <list>          <list>            <list>
## GPGQECA                      0.505943       0.0236133         0.0466718
## CAVCGWALPHNRM                0.819962       0.0111598         0.0136102
## CAVCGWALPHNRMQAL              3.18841        0.130744         0.0410059
## CAVCGWALPHNRMQALTSCE          1.51203        0.598316          0.395703
## AVCGWALPHNRM                 0.138914      0.00401256         0.0288852
## ...                               ...             ...               ...
## ATERYLHVRPQPLAGEDPPAYQARL    0.882916       0.0634907         0.0719102
## ERYLHVRPQPLAGEDPPAYQ         0.841277       0.0671491         0.0798181
## ERYLHVRPQPLAGEDPPAYQARL       2.11026        0.172771          0.081872
## RYLHVRPQPLAGEDPPAYQARL        1.65969        0.120275         0.0724686
## LQKLTEEVPLGQSIPRRRK           4.48223        0.196504         0.0438406
##                               pvals       fdr ebayes.pvals ebayes.fdr
##                           <numeric> <numeric>    <numeric>  <numeric>
## GPGQECA                    0.516181  0.652643     0.479007   0.612683
## CAVCGWALPHNRM              0.416403  0.558589     0.406959   0.545920
## CAVCGWALPHNRMQAL           0.148709  0.320744     0.123017   0.300707
## CAVCGWALPHNRMQALTSCE       0.286210  0.499732     0.237829   0.441045
## AVCGWALPHNRM               0.728272  0.801099     0.708845   0.779730
## ...                             ...       ...          ...        ...
## ATERYLHVRPQPLAGEDPPAYQARL  0.400605  0.557804    0.3563859   0.502596
## ERYLHVRPQPLAGEDPPAYQ       0.410930  0.558053    0.3658748   0.509446
## ERYLHVRPQPLAGEDPPAYQARL    0.219967  0.443454    0.1817334   0.386882
## RYLHVRPQPLAGEDPPAYQARL     0.267115  0.481684    0.2263725   0.436859
## LQKLTEEVPLGQSIPRRRK        0.101671  0.266280    0.0814249   0.246599
##                           empirical.fdr fitcomplete
##                               <numeric>   <integer>
## GPGQECA                        0.795215           1
## CAVCGWALPHNRM                  0.778284           2
## CAVCGWALPHNRMQAL               0.778284           3
## CAVCGWALPHNRMQALTSCE           0.778284           4
## AVCGWALPHNRM                   0.816964           5
## ...                                 ...         ...
## ATERYLHVRPQPLAGEDPPAYQARL      0.778284         106
## ERYLHVRPQPLAGEDPPAYQ           0.778284         107
## ERYLHVRPQPLAGEDPPAYQARL        0.778284         108
## RYLHVRPQPLAGEDPPAYQARL         0.778284         109
## LQKLTEEVPLGQSIPRRRK            0.778284         110
which(funresdAb25_1@results$ebayes.fdr < 0.05)
##           IQLRESLEPDA RESLEPDAYALFHKKLTEGVL         YALFHKKLTEGVL 
##                    36                    42                    43 
##       REQLEATCPQCHQTF           EATCPQCHQTF       MYLQENGIDCPKCKF 
##                    52                    53                    65 
##     YLQENGIDCPKCKFSYA      LQENGIDCPKCKFSYA 
##                    68                    70

We can plot these kinetics to see what is happening. This allows us to visualise region of protection and deprotection, potentially identifiying the epitope.

(res_all@statmodels[[36]]@vis + 
res_all@statmodels[[42]]@vis  + 
res_all@statmodels[[43]]@vis  + 
res_all@statmodels[[65]]@vis  + 
res_all@statmodels[[68]]@vis  + 
res_all@statmodels[[70]]@vis  + 
res_all@statmodels[[52]]@vis  + 
res_all@statmodels[[53]]@vis ) + plot_layout(guides = 'collect')

We can make a Manhattan plot to better specially visualise what’s happening.

#We need to provide an indication of "difference" so we can examine deprotected
# or prected regions
diffdata <- assay(newqDF)[,6] - assay(newqDF)[,3]


sigplots <- manhattanplot(params = funresdAb25_1,
                          sequences = HOIP$Sequence, 
                          region = HOIP[, c("Start", "End")],
                          difference = diffdata,
                          nrow = 1)
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
sigplots[[1]] + plot_layout(guides = 'collect')

We can visualise this in a peptide plot which helps us understand the nature of the overlap

fpath <- system.file("extdata", "HOIP.txt", package = "hdxstats", mustWork = TRUE)
HOIPfasta <- readAAStringSet(filepath = fpath, "fasta")

scores <- funresdAb25_1@results$ebayes.fdr

out <- plotEpitopeMap(AAString = HOIPfasta[[1]],
                      peptideSeqs = unique(HOIP$Sequence),
                      numlines = 2,
                      maxmismatch = 1,
                      by = 1,
                      scores = 1 * (-log10(scores[unique(HOIP$Sequence)])  > -log10(0.05)) + 0.0001,
                      name = "significant")
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
out[[1]]/(out[[2]]) + plot_layout(guides = 'collect') & theme(legend.position = "right")

We can further visualise this a barcode of particular residues, here we use residue level averaging to obtain results at the residue level.

scores <- funresdAb25_1@results$ebayes.fdr
out2 <- plotEpitopeMapResidue(AAString = HOIPfasta[[1]],
                              peptideSeqs = unique(HOIP$Sequence),
                              numlines = 2,
                              maxmismatch = 1,
                              by = 5,
                              scores = scores[unique(HOIP$Sequence)],
                              name = "-log10 p value")
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
out2[[1]]/out2[[2]]  + plot_layout(guides = 'collect') & theme(legend.position = "right")

We can also plot multiple residue maps on the same plot so that we can compare different antibodies.

scores <- funresdAb25_1@results$ebayes.fdr
avMap25_1 <- ComputeAverageMap(AAString = HOIPfasta[[1]],
                               peptideSeqs = unique(HOIP$Sequence),
                               numlines = 2, maxmismatch = 1,
                               by = 10, scores = scores[unique(HOIP$Sequence)],
                               name = "-log10 p value")

## generate results from other dAB
cn <- new.colnames[c(1:3,19:21)]
HOIP_wide_nrm <- data.frame(HOIP_wide)
HOIP_wide_nrm[,c(2:4)] <- HOIP_wide_nrm[,c(2:4)] - HOIP_wide_nrm[,c(2)]
HOIP_wide_nrm[,c(20:22)] <- HOIP_wide_nrm[,c(20:22)] - HOIP_wide_nrm[,c(20)] 

newqDF2 <- parseDeutData(object = DataFrame(HOIP_wide_nrm),
                        design = cn,
                        quantcol = c(2:4,20:22),
                        rownames = HOIP_wide$Sequence)

res_all2 <- fitUptakeKinetics(object = newqDF2[,1:6],
                             feature = rownames(newqDF2[,1:6])[[1]],
                             start = list(a = NULL),
                             formula = value ~ a * (1 - exp(-0.07*(timepoint))), maxAttempts = 1)
## Warning in max(data$value): no non-missing arguments to max; returning -Inf
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
funresdAb27_2 <- processFunctional(object = newqDF[,1:6],
                                   params = res_all2)
## Warning in stats::optim(x = c(9.29818825569032e-05, 0.090283790910004, 0.166406795432091, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in densfun(x, parm, ...): NaNs produced
## Warning in stats::optim(x = c(0.00297145686513731, 0.032686171739902, 0.0181780887330962, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in densfun(x, parm, ...): NaNs produced
scores <- funresdAb27_2@results$ebayes.fdr
# compute average map
avMap27_2  <- ComputeAverageMap(AAString = HOIPfasta[[1]],
                                peptideSeqs = unique(HOIP$Sequence),
                                numlines = 2,
                                maxmismatch = 1,
                                by = 10,
                                scores = scores[unique(HOIP$Sequence)],
                                name = "-log10 p value")

# set rownames
rownames(avMap25_1) <- "dAb25_1"
rownames(avMap27_2) <- "dAb27_2"

# store in a list
avMap <- list(avMap27_2 = avMap27_2,
              avMap25_1 = avMap25_1)

#plotting
out3 <- plotAverageMaps(avMap, by = 20)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
out3[[1]]/out3[[2]]  + plot_layout(guides = 'collect') & theme(legend.position = "right")

9 Protection/deprotection maps

Users may also wish to colour according to protection or deprotections of hdx. The following functions determine whether there is protection or not by computing deuterium difference. The users pass the argument for which to compute the differences. The user pass can any scores though it is useful to look at the FDR provided by previous analysis. As is typical blue indicates protection whilst red deprotection. The binding epitope is very clear from these diagrams.

scores <- funresdAb27_2@results$ebayes.fdr
hdxdiff1 <- hdxdifference(object = newqDF2,
              AAString = HOIPfasta[[1]],
              peptideSeqs = unique(HOIP$Sequence),
              numlines = 2,
              maxmismatch = 1,
              by = 10,
              scores = scores[unique(HOIP$Sequence)],
              cols = c(2,5), # columns of object to use for differences
              name = "-log10 p value (signed)")
## Warning in hdxdifference(object = newqDF2, AAString = HOIPfasta[[1]],
## peptideSeqs = unique(HOIP$Sequence), : NaNs produced
scores <- funresdAb25_1@results$ebayes.fdr
hdxdiff2 <- hdxdifference(object = newqDF,
              AAString = HOIPfasta[[1]],
              peptideSeqs = unique(HOIP$Sequence),
              numlines = 2,
              maxmismatch = 1,
              by = 10,
              scores = scores[unique(HOIP$Sequence)],
              cols = c(2,5),
              name = "-log10 p value (signed)")
## Warning in hdxdifference(object = newqDF, AAString = HOIPfasta[[1]], peptideSeqs
## = unique(HOIP$Sequence), : NaNs produced
dMap <- list(hdxdiff1$diffMap, hdxdiff2$diffMap)
out4 <- hdxheatmap(averageMaps = avMap, diffMaps = dMap)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
out4[[1]]/out4[[2]]  + plot_layout(guides = 'collect') & theme(legend.position = "right")

10 Exporting results

We suggest saving files in .rds format. However, if you wish to export data to .csv format this can be done as in the following exampe.

# not run
write.csv(funresdAb25_1@results, file = "myhdxresuts.csv")